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Abstract 

The microscopic dynamics of one-dimensional self-gravitating many-body 
systems is studied. We examine two courses of the evolution which has the 
isothermal and stationary water-bag distribution as initial conditions. We in- 
vestigate the evolution of the systems toward thermal equilibrium. It is found 
that when the number of degrees of freedom of the system is increased, the 
water-bag distribution becomes a quasi-equilibrium, and also the stochastic- 
ity of the system reduces. This results suggest that the phase space of the 
system is effectively not ergodic and the system with large degreees of freedom 

approaches to the near-integrable one. 
05.45.+b, 98.10.+Z, 03.20.+i, 95.10.Ce 
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I. INTRODUCTION 



A self-gravitating many body system, which contains many particles interacting with mu- 
tual gravity, is an idealized model of wide classes of astronomical objects, such as globular 
clusters, elliptical galaxies, and clusters of galaxy. Dynamics of the system has two char- 
acteristic aspects. One is the microscopic dynamics, which is concerned with the motions 
of the individual particles. The other is the macroscopic dynamics, which deals with the 
averaged quantities. Evolution of the real system is determined exactly by the microscopic 
dynamics. However, the number of the particles contained in the system is so large, e.g. 
10 5 for globular clusters and 10 11 for elliptical galaxies, that we can treat only macroscopic 
quantities by the statistical method, in practice. 

In many cases which are successfully treated by statistical mechanics, which include gas 
and liquid systems, the particles interact with close neighbours, then the local equilibrium 
is determined by the state of the neighbours. Thus if the macroscopic quantities are defined 
by averaging over the scale much larger than the range of the force, their evolution is 
independent of the microscopic dynamics and is governed only by themselves . In these 
systems the statistical mechanics is easily applicable. 

On the other hand, in the self-gravitating systems, the motions of the individual particles 
are governed by the summation of the forces from all the other, because the gravity is 
long range force. Thus the macroscopic dynamics is not decoupled with the microscopic 
dynamics. Therefore we must study how the microscopic dynamics influences the evolution 
of the macroscopic quantities in the self-gravitating systems. 

The microscopic dynamics is described as a trajectory in the T-space (iV-body phase 
space). If the system has some integrals of motion, the trajectory is confined in the subspace 
which conserves the integrals. The system is defined as ergodic if the time average of any 
dynamical quantities are equal to the spatial average over the subspace. This subspace is 
referred to as the ergodic region. In this case the time average does not depend on the initial 
condition of the system but only on the integrals, and thus it is time independent. This 



2 



gives the statistical or micro-canonical equilibrium. At the equilibrium, energy is equally 
distributed to all degrees of freedom, which is called equipartition. 

The mixing system, which is included in the ergodic system, has the property that a 
small but finite part of the phase space volume spreads over whole ergodic region by means 
of coarse-graining. This property is closely related to relaxation of the system, because 
information of the initial state is lost and the micro-canonical equilibrium is realised. Any 
time-correlation disappears after the trajectory is mixed over the ergodic region [for more de- 
tailed reviews, see, e.g. We refer to the time that mixing is realized as the microscopic 
relaxation time. 

Though ergodicity is the fundamental basis of statistical mechanics, not all systems are 
ergodic. For example, the phase space of a near-integrable Hamiltonian system contains 
both stochastic region and regular region (tori). A trajectory initially located on a torus 
stays on it forever. This system is not ergodic over the subspace where the total energy of 
the system is constant. Furthermore, it is shown that there exist stagnant layers around 
tori and a trajectory in the stagnant layer stays there for a long time ||. In this system, the 
trajectory in the stagnant layer moves in the stochastic region after a long interval of time, 
and then they could enter another stagnant layer. Thus, in the system, transient stationary 
states can emerge in some era. We refer to them as quasi- equilibria. 

Lynden-Bell applied the concept of ergodicity to the evolution of elliptical galaxies, and 
derived the unique equilibrium state ||. However, observations of elliptical galaxies and 
numerical simulations show disagreement with the Lynden-Bell distribution 0-0. In par- 
ticular, elliptical galaxies are believed to be triaxial in the shape and anisotropic in the 
velocity dispersion. This stationary state seems to suggest the existence of additional inte- 
grals, which conserves the anisotropy, though it is unclear whether the state is a stationary 
state induced by the additional integrals or only a transient state approaching to the sta- 
tistical equilibrium. These facts suggest that elliptical galaxies are not ergodic while it is 
generally believed that the self-gravitating many body systems are chaotic and so ergodic. 
So it is very important and interesting to study ergodicity of the self-gravitating many body 
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systems in order to analyse the present dynamical structures of elliptical galaxies, and more- 
over to examine applicability of usual statistical mechanics to the self-gravitating many-body 
systems. 

To study their ergodicity, we employ one-dimensional system. Because in one-dimension 
the phase space is compact, which makes the system tractable in considering ergodicity. 
Another reason is that the force law is very simple and thus the evolution of the system can 
be followed numerically with a good accuracy by using exact code ||. Though the force law 
in one-dimensional system is different from that in the three-dimensional system, we can 
study the properties induced by long range forces even in the one-dimensional systems. 

Hohl |9||| firstly suggested that the one- dimensional system consisting of N identical 
plane-parallel sheets should relax on the order of N 2 t c , where t c is the characteristic time 
which is approximately the time for a member to traverse the system. However, Wright, 
Miller, and Stein |H| asserted that some initial states does not approach to micro-canonical 
distribution after 2N 2 t c . Thus the Hole's conjecture was questioned. In succeeding pa- 
per |Ti|-|r7|], it was shown that the evolution of the system greatly depends on the initial 
condition. Some initial state appears to relax on the time scale of Nt c , which is much 
shorter than Hohl's prediction ]T1|,T^,I^]. Then the complicated features of relaxation in the 



self-gravitating systems were recognized. 

Severne and Luwel |19| suggested that there are three phase in relaxation. If the initial 



state is far from equilibrium, violent oscillation of mean field gives rise to the violent relax- 
ation j| for the first several oscillation. After the system almost virialised, remaining small 
fluctuation of gravitational field causes the change of the individual particle energies. They 
called this era as the collisionless mixing phase. In astrophysics the evolution in the violent 
relaxation and collisionless mixing phase are often referred to dynamical evolution. After 
that, the collisional relaxation phase takes place, in which the particle interactions tend to 
drive the system towards the microscopic thermal equilibrium. This is the thermal evolu- 
tion. Luwel and Severne |TB| showed that the collisionless mixing occurs in the stationary 



water-bag distribution. The distribution is one of the stationary solutions of the collisionless 



Boltzmann-Poisson equations. Thus it does not change its macroscopic distribution in colli- 
sionless phase. Luwel and Severne did not study the evolution in the collisional phase. If the 
system is ergodic, it is expected that the stationary water-bag distribution is transformed 
into the microcanonical distribution after the collisional relaxation becomes effective. 

While the above studies are concerned with the evolution of the macroscopic quantities, 
the microscopic dynamics, especially the ergodicity of the system also has been studied by 
some authors. For small N (N < 10) it is shown that the systems are "ergodic" for N > 5 



22|. For 11 < iV < 20, Reidl and Miller |T^] found that two different periodic orbits 
have the different maximum Lyapunov exponents. From this fact they asserted that the 
trajectories of the two orbits in the phase space covers different regions, which are separated 
from each other. 

However, this conclusion is doubtful because the measure of the trajectory of the periodic 
orbit in the phase space is negligible then the study gave no information of the whole 
phase space dynamics. Further, the numbers of the sheets are too small to speculate the 
macroscopic relaxation in larger N. Thus we performed numerical simulations for N = 32, 
128, 512, and examine the dependence of the dynamics on the number of the sheets (i.e., 
the degrees of freedom of the system). 

In order to get more information of microscopic dynamics, we have investigated the 
time correlation of the fluctuation of the individual particle energies, the equipartition of 
the individual particle energies and the convergence of the maximum Lyapunov exponent. 
In general, the motion in an ergodic region shows fast decay of correlation. On the other 
hand, in the near integrable Hamilton system an essential feature of the diffusion process 
(including the Arnold diffusion) is the appearance of long time tails of the correlation or 
the enhancement of the diffusion mode with the zero frequency, e.g., the power spectrum 
density (PSD) function S(f) satisfies, 

s(f)~r v , (/«i) (i) 

where / stands for the frequency and u a positive constant (2 > v > 1) |3|j23||. Therefore the 



time correlation or the PSD is a good tool to understand the ergodicity of the phase space, 
besides the equipartition of the energies and the convergence of the maximum Lyapunov 
exponents. 

In section [TT] we describe the model and the initial conditions. Three quantities which 
we employ to analyse the system are explained in section [ITJ and results of the numerical 
simulations are given in section |IV|. We devote section |V| to the conclusions and discussions. 



II. DESCRIPTION OF THE MODEL 

The one-dimensional self-gravitating system consists of N identical mass sheets, each of 
uniform mass density and infinite in extent in the (y, z)-plane. We call the sheet as particle 
in this paper. The particles are free to move along and accelerate result of their 

mutual gravitational attraction. The Hamiltonian of the system is given by 

# =yX>? + (27rGm 2 ) J>,. - Xi \, (2) 

i=l i<j 

where m, V{, and Xi are the mass (surface density), velocity, and position of the ith par- 
ticle, respectively. Since the gravitational field is uniform, the individual particles moves 
parabolically, until they intersect with the neighbours. Thus the evolution of the system can 
be followed by connecting the parabolic motions. When an encounter occurs between two 
particles, they pass freely through each other. 



Our code is the application of that for one-dimensional sheet plasma [24| to the self- 



gravitating systems, which is very similar to the "new exact code" referred in ref | 19| , but 
the exchange of particles at an encounter is arranged by the heap sort algorithm. During 
the integration, time is measured in the unit 

t c = (47rG Pav )~ 1/2 , (3) 

where p av is the mass divided by the width of the distribution at the initial time. 

All calculation were performed in double precision (16 significant figures) on DEC station 
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3000AXP and a SONY NEWS 5000. The total energy was conserved to better than one 
part in 10 13 . 

We examined two courses of evolution, IT and WB, which begin with the isothermal and 
water-bag distributions as their initial condition, respectively. In the water-bag distribution, 
all particles are randomly distributed in a rectangle region in the /i-space to form a uniform 
averaged phase density (Fig.|l]). This is not a exact stationary solution of the collisionless 
Boltzmann-Poisson equations because the shape is a rectangle, but the difference is very 
small and then after several oscillations it is transformed into the stationary water-bag 
configuration. The isothermal distribution is given by 

0(77) = 7r~ 1/2 exp(-?7 2 ) (velocity), (4) 

P(0 = ^ secn2 ^ (position), (5) 

where 

V = (v/2)(3M/E)^ 2 , (6) 

and 

f = (3ttGM 2 /2E)x, (7) 

M and E represent the total system mass and total energy, respectively(Fig. |2|). It is a 
stationary solution of the collisionless Boltzmann-Poisson equations. In all our simulations, 
the total mass and energy set to 1 and 1/4, respectively. The typical period of oscillation of 
a particle is 2nt c . 

We choose the initial distributions which are dynamically stable, because our studies are 
concentrated on the thermal evolution (collisional relaxation) of the system. If the system 
is ergodic, IT and WB should coincide after the relaxation time. Therefore we compare the 
behaviours of IT and WB for a long time. 

III. ANALYSES 
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A. Equipartition 



The specific energy (energy per unit mass) £*(£) of ith particle is given by 

1 N 
= ~ v Kt) + 2irGmY,\xj(t) ~ (8) 

^ 3=1 

If the evolution of the system is ergodic in the T-space, the long time average of the specific 
energy takes a unique value for all i, i.e. 

r t =]^^J%i(t)dt = eo = 5E/3. (9) 

The degree of the deviation from the equipartition is measured by the quantity, 



A(t)= £ -\ 



1 N 

m^E(^)-^o) 2 , (io) 



where £l(t) is the averaged value until t. In the numerical scheme, £j(t)s are sampled at 
every At = 0.78125t c , and the average is defined simply by the summation of the samples 
divided by the number of the samples. 

If the system is mixing, we can estimate the way of the evolution of A(t). In the system, 
time correlation disappears in a finite time (the relaxation time). In this time, a trajectory 
in the T-space visits almost every point in the ergodic region, thus the A(t) almost vanishes 
but gives a small fluctuation. Since the correlation vanishes after the relaxation time, A(t) 
decreases as t~ 1//2 for the time longer than the relaxation time, according to the central limit 
theorem. Therefore the evolution of A(t) is one of good tests of mixing property of the 
system. 

The initial value of the deviation from the equipartition, A(0), depends on the initial 
distribution. Figure [3] shows the cumulative distribution of the specific energy, 0(e), which 
is defined by 

u(e) = (l/N)N(ei < e), (11) 

where N(Ei < e) is the number of particles with £j < e. The isothermal and water-bag 
distribution are represented by the solid and dashed lines, respectively. Also the energy 
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at the equipartition, e , is shown by dotted line. Since the isothermal distribution spreads 
wider in the energy space than the water-bag, A(0) of the isothermal distribution is larger 
than that of the water-bag distribution. 



B. Power spectrum density 

The power spectrum density (PSD) S(f) is given by 



N 

S(f) = (l/N)Y,\C>(f)\ 2 , (12) 



i=i 

where Ci(f) is Fourier transform of £j(i), i.e. 



*(f) = / aUle'^df. (13) 



The PSD, on the other hand, is the Fourier transform of the autocorrelation function, thus 
the long time correlation gives rise to a peak at small /. In order to obtain PSD numerically, 
a sequence of "locally averaged" energy {(e) (ti = T /n), (e)(ia — 2T /n), . . . , (e)(t n = T )} 
are sampled, where 

(e)i(tj) = - f 3 Ei{t)dt. (14) 
ti Jtj-i 

The maximum integration time To determines the minimum frequency / m j n = 1/Tq, and 
the interval of sampling determines the maximum frequency / max . This local averaging 
suppresses the higher frequency modes than / max , thus it prevents the "aliasing", which 



means that the higher frequency modes fall into the interval of / < / max ||25|| . The procedure 
of the local averaging, eq. (]T4|), is the same as that described in the previous section. We 
integrated the motion to T = 10 7 t c for iV =32 and 128, hence /m in = lO -7 ^" 1 . The 
number of samples n is limited to 256 due to the computer's ability, hence / max = 2.56 x 
10 _5 t~ 1 . In order to obtain higher /, we gathered n samples with shorter time interval, 
{(e<)(T /10n), (ei){T /10n), (e i )(2r /10n), ■ ■ • , (e<)(r /10)}, which covers the range of the 
frequency from lO -6 ^" 1 to 2.56 x 10 — 4 t ( 7 1 . If the amplitude of oscillation with the frequency 
KT 6 ^ 1 < f < 2.56 x KT 5 ^ 1 for t > 10% is the same as during t < 10%, two S(f) merge 
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into one curve. In this way we extended the range of the frequency beyond 10 — 2 t c 7 1 . The 
system with N = 128 is the exception that the values of S(f) of T = 10 7 t c and T < 10 6 t c , 
for 10 6 t^T 1 < / < 2.56 x 10 _5 t~ 1 are different (Fig. [7]). This is discussed in the next section. 
Since the PSD is the Fourier transform of the autocorrelation function 

= je(t + t')e(t')dt', (15) 

we can find the time scale on which the correlation decays. For example, the Brownian 
motion, u(t), subjected by the Langevin equation, 

^ + lu = R(t), (16) 

where 7 is the friction coefficient and R(t) is a random force (white noise), gives the PSD 
of the Lorentz distribution, 

S (f) « J^-l- ( 17 ) 

For / ^> 7, S(f) oc f~ 2 , which describes the short time scale nature of the system. In 
this phase, the random force R(t) is dominated and the dispersion ((u(t) —u(0)) 2 ) increases 
proportional to t. In other words, it is the diffusion phase from its initial value. On the 
other hand, for /C7, S(f) is constant. Thus any correlation disappears in the time scale 
7 _1 , and then the fluctuation becomes like a thermal noise. Since this time scale is closely 
related to the relaxation time, the PSD gives us a complemental information about the 
mixing property of the system. 

We should mention more about the PSD with a peak or divergence as / goes to zero, 
e.g. S(f) oc which means the existence of a long time correlation. Such a PSD is 
often observed in a near-integrable system. In the system, the phase space is divided into 
stochastic region and regular region (tori). It is well-known as the Arnold diffusion that it 
takes a long time for a trajectory to travel across the web of the tori. This slow diffusion 
gives rise to a long time correlation |3|]. 
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C. Lyapunov exponent 



The evolution of the system can be described a trajectory in the T-space. The Lyapunov 
exponents of a given trajectory characterized the mean exponential rate of divergence of 
trajectories surrounding it [a general discussion of the Lyapunov exponents is provided by 
Lichtenberg and Lieberman H], 

The maximum Lyapunov exponent is defined as 

X = lim -In^r, (18) 
*(o)-ot d(0Y V ; 

where d(t) and d(0) are the separations in the T-space between two nearby orbits at times t 



and 0, respectively. The numerical procedure follows mostly Shimada and Nagashima |26 
where d(t) is determined by the linearized equations of motion: 

Aii = Avi, (19) 

N 

Avi = -A7cGm^2(Axi - Axj)S(xi - xj), (20) 



d 



N 

J2(Ax* + Av?), (21) 



where Axi and Avi are the first order deviation of the position and velocity. Since d(t) 
diverges exponentially, when the separation d becomes 10 5 times the initial value, the sepa- 
ration is rescaled to d(0). 

The Lyapunov exponent is originally defined by the limit of infinite t, but here we use 
"time-dependent Lyapunov exponent" defined by 

w = i h H (22) 

If the system is mixing, the time average becomes roughly the same as the space average over 
ergodic region, in the relaxation time. Thus the convergence of X(t) gives another measure 
of relaxation. 
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IV. RESULTS 
A. iV=32 

Figure |] shows the time evolution of the deviation from the equipartition, A(£) (bold 
lines) and the Lyapunov exponents, \(t) (light lines). The PSD, S(f), are shown in Fig. |5|. 
The curves of IT and WB are represented by solid and dashed lines, respectively. 

The features of the curves are not so simple as we expected in the previous section. A(t) 
begins to decrease at t ~ 3 x 10 3 t c , but does not show monotonic decrease like t^ 1 ^ 2 . From 
t ~ 3 x 10 3 t c , A(£) of IT and WB coincide, but at t ~ 10 4 t c , IT begins to increase again 
and the difference between IT and WB lasts until t ~ 3 x 10 6 t c . After that, though the 
range of t is not so long, A(i) seems to decrease as t~ x l 2 . At this time, t ~ 3 x 10 6 t c , also 
the convergence of \{t) becomes quite well. From the PSD (Fig. |5]), it is found that S(f) 
is constant for / <C 10 — 5 t ( 7 1 , which gives the time scale of disappearance of correlation of 
t ~ 10 5 t c . Therefore it is probable that the system relaxes on this time, t ~ 10 6 t c . 

From Fig. f|, there seems to exist another time scale, t ~ 10 3 t c , at which A(t) begins to 
decrease. For shorter period, t 10 3 t c , the PSD indicates that the variation of £i(t) is well- 
approximated by the random walk because S(f) oc f~ 2 . For t > 3 x 10 3 t c , the trajectory 
in the T-space travels farther so that the individual particle energies vary widely over the 
energy space. However, the long time correlation seems to last for several 10 4 t c , because 
A(£) does not decrease monotonically like t~ x l 2 and convergence of A(t) is not good. This 
correlation is not swept away by the diffusion due to the random walk. 

From these facts, it is speculated that the system has a part of mixing property with 
unique ergodic region, but the phase space is not uniform, and contains many structures, 
maybe a set of ruined tori, which yield the long time correlation for t < 10 5 t c . 
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B. jV=128 



Figure || and are the same as Fig. | and || but for N = 128. The PSD for WB is 
that S(f) is separated into two part, lO -7 *" 1 < f < 2.56 x lO -5 *" 1 and lO" 6 ^ 1 < f < 
2.56 x lO -2 ^" 1 . The former is calculated from the sample of To = 10 7 t c and n = 256 (see 
I TTD3) . The other is calculated from that of T < 10 6 t c , thus this part does not contain the 



behaviour during 10 6 t c < t. The difference between the two means that evolution with a 
period 10 6 t c later than t = 10 6 t c is different from that of < t < 10 6 t c . The transition of 
state at t ~ 10 6 t c is found in Fig. |5|. A(t) of WB decreases as t~ 1 / 2 until t < 10H C and then 
suddenly increases to approach to IT. Also X(t) shows the transition at t ~ 2 x 10 6 t c . After 
that A(t) and X(t) of WB show tendency to approach to that of IT. Therefore the system 
is expected to relax on the time scale of 10 7 t c . 

Another interesting feature of this system is the behaviour of WB for t < 10 6 t c . A(t) 
shows a good agreement with t _1//2 from t ~ 10 2 t c , and X(t) converges to the local value, 
5.5 x 10~ 2 , while that of IT converges to 6.2 x 10~ 2 . Furthermore, S(f) calculated from the 
sample of T < 10 6 t c shows almost flat spectrum for / < 10 _3 t c 7 1 . These features indicate 
that WB relaxes in t ~ 10 3 t c to a quasi-equilibrium which lasts until t ~ 10 6 t c . Figure |8| 
shows the cumulative energy distribution, ^(e), of IT at the beginning, t = 0, (solid line), 
WB at the beginning (dashed line), and WB at t = 10 6 t c (dotted line). The distribution 
of WB at t = 10 6 t c is very close to the initial distribution compared with the isothermal 
distribution. Therefore the quasi-equilibrium has, in fact, the water-bag distribution. On 
the other hand, the difference between IT and WB disappears after the transition. Figure 
P shows v(e) of IT (solid line) and WB (dashed line) at t = 10 7 t c . 

The evolution of IT is similar to that of iV = 32. There is a long time correlation, which 
is found from the evolution of A(t) and A(t). Especially for iV = 128, PSD shows clear bend 
at / = 2 x lO" 4 ^ 1 , and the power of / for lO" 6 ^ 1 < f < lO^H" 1 is less steep than f~ 2 . 
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C. N = 512 



Figure 10 and 11 show the same as Fig. |] and [5] but for N = 512, respectively. 

The PSD indicates that the behaviour in a shorter time scale than 10 3 t c is the same for 
IT and WB. It is the diffusion process of the individual particle energy from its initial value. 

Since the integration time is limited to 10 6 t c from the constraint of computer's ability, 
the transition of WB from the water-bag to the isothermal distribution, which is expected 
to occur at later than 10 6 t c , is not observed. The behaviour of WB for t < 10 6 t c is quite 
similar to that for N = 128. Thus water-bag distribution is one of quasi-equilibrium and 
the relaxation time to the quasi-equilibrium is about several 10 3 t c . 

Increasing the number of the particles for IT strengthens the property that disappearance 
of time-correlation which is found by the bend of the PSD occurs the later for the larger N. 
In fact, the frequency at which the feature of PSD changes from flat to inclined is lO -5 ^" 1 
for N = 32 and 10 _6 t ( 7 1 for iV = 128. Together with the fact that A(t) does not decrease 
as t -1 / 2 , we found that the relaxation time, at which the trajectory covers almost whole 
ergodic region, becomes longer for larger N. 



V. CONCLUSIONS AND DISCUSSIONS 

We performed iV-body simulations of the evolution of one-dimensional self-gravitating 
systems for different N. We examined two courses of evolution, IT and WB, which have the 
initial distributions of the isothermal and water-bag, respectively. We obtained the following 
lines of conclusion. 

1. These systems relax to the isothermal distribution. The time at which the difference 
between IT and WB disappears is about 10 6 t c for N = 32, 10 6 t c for N = 128, and 
longer than 10 6 t c for N = 512. 

2. As N increases, the water-bag distribution becomes a quasi-equilibrium, which yields 
ergodic properties apparently. The microscopic relaxation time of the water-bag, at 
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which the time correlation disappears, is 10 3 t c for N = 128 and 10 4 t c for iV = 512. 
The water-bag distribution lasts stably for a period orders of magnitude longer time 
than the time that the correlation disappears. 

3. Since the isothermal distribution is the micro-canonical equilibrium, IT shows no evo- 
lution macroscopically. As one can see in Figures of PSD, the microscopic relaxation 
time is, however, much longer than that of the water-bag quasi-equilibrium. For 
N = 32, the microscopic relaxation time is about 10 5 t c and 10 6 t c for N = 128. For 
iV = 512, we did not observe the microscopic relaxation in 10 6 t c . The transition of the 
water-bag to the isothermal distribution occurs in the same time scale. 

4. Convergence of the maximum time-dependent Lyapunov exponents to the different 
values for IT and WB indicates that the phase space is separated into at least two 
ergodic-like regions which give the quasi-equilibria. 

In the limit of iV — > oo the evolution of the system, both the isothermal and water-bag 
distribution are the stationary solution of the collisionless Boltzmann-Poisson equations. 
The distribution function which depends only on the individual particle energy is always 
the stationary solution of the collisionless Boltzmann equation. Among these solutions 
many satisfy the Poisson equation, so many class of analytic stationary solutions for the 
collisionless Boltzmann-Poisson equations exist. These solutions are expected to be also 
quasi-equilibria, and they should depend on the initial condition. In fact, we have examined 
another initial condition, in which all particles have the same energy. The distribution of 



particles is shown in Fig. |12| Since this distribution is unstable, the ellipse in the phase 
space is broken rapidly, then the system approaches not to the isothermal distribution but 
to the water-bag distribution. So in this case the final results are similar to the one for 
WB mentioned before. Furthermore we are presently studying the dependence of the initial 
condition in more details. The result will be described in the succeeding paper. 

From these facts, it is obvious that the phase space is not homogeneous but contains many 
substructures which induce the quasi-equilibria. Though the origin of the substructures is 
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not clear, we have a speculation that tori play an important role. This idea conflicts with 
the usual idea that chaos is strengthen as the number of degrees of freedom increases, thus 
the system approaches to ergodic. The latter idea is supported by some numerical results, in 
which the volume of the KAM tori is reduced as the number of degrees of freedom increases 
27|,2B|. However, we believe our idea, because of the following reasons: in the limit of 



N — > oo, the distribution function which depends only on the individual particle energy is 
alway the stationary solution of the collisionless Boltzmann equation. These solutions form 
tori because all phase elements in the phase space (i.e. all degrees of freedom) conserves their 
energy, and the motion which has the same number of integrals of motion as the degree of 
freedom is a torus. Thus in the limit the system contains many tori, which seem to produce 
the multi-ergodicity ||. The other reason is that, in fact, chaos is getting weaker as N 
increases. Figure [13| shows the converged value of the maximum time-dependent Lyapunov 
exponents of IT and WB for various N. That of IT decreases as iV~ 1/5 and of WB as iV" 1 / 4 . 
This fact suggests that the system approaches to near-integrable system as N increases. 

The property that chaos becomes weaker as the number of degrees of freedom of the 
system increases seems due to the long range force, because the force acting on a particle 
is the summation of the force of all the others, thus the effect of close encounter, which 
is usually believed to be the origin of chaos, becomes weaker as the system population 
increases. The results for N < 10 that the stochasticity of the systems are strengthened 
|20| p2[ seem to be due to the small number of particles in the system. In such systems, the 
properties of the long range force are not prominent. Thus it is inferred that there exists 
the number of degrees of freedom which maximizes the stochasticity of the system between 
10 and 30. 

In the three-dimensional system, the properties described above are expected to hold, 
because the force is also long range. In fact, elliptical galaxies have the stationary state 
which is triaxial in shape and anisotropic in the velocity dispersion. Also the stationary 
states depend on the initial conditions []51||. This state is completely different from what 
Lynden-Bell predicted by using ergodic hypotheses [§J. This state can be explained as the 
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quasi-equilibrium mentioned in this paper. 
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FIGURES 



FIG. 1. The water-bag initial distribution in /i-space. 

FIG. 2. The isothermal initial distribution in jU-space. 

FIG. 3. The cumulative specific energy distribution of the isothermal (solid line) and water-bag 
(dashed line) distribution for N = 512. The dotted line indicates the energy at equipartition. 

FIG. 4. Evolution of the deviation from the equipartition, A(t) (bold lines), and the maximum 
Lyapunov exponent, X(t) (light lines) for TV = 32. The solid and dashed lines represents IT and 
WB, respectively. The unit of time is t c . 

FIG. 5. Power spectrum densities of IT (solid line) and WB (dashed line) for N = 32. The 
unit of frequency is i" 1 . 

FIG. 6. Same as Fig. 4, but for N = 128. 

FIG. 7. Same as Fig. 5, but for N = 128. 

FIG. 8. Cumulative energy distributions, 0(e), of IT at the beginning, t = 0, (solid line), WB 
at the beginning (dashed line), and WB at t = W e t c (dotted line). 

FIG. 9. Cumulative energy distributions, 0(e), of IT (solid line) and WB (dashed line) at 
t = W 7 t c . 

FIG. 10. Same as Fig. 4, but for N = 512. 
FIG. 11. Same as Fig. 5, but for N = 512. 
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FIG. 12. Distribution of particle in a stationary solution of collisionless Boltzmann-Poisson 
equations. All the particles have the same energy. 

FIG. 13. The maximum Lyapunov exponents of IT (solid line) and WB (dashed line) after 
convergence. 
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